
/*
 *  interface.cpp
 *  EPPI-0.0
 *
 *  Created by Juan Carlos Garcia Orden on 15/11/07.
 *
 */

#include "interface.h"

Interface::Interface()
:	lProc(0),
	rProc(0),
	localProc(),
	relatedInterfaceNum(0),
	isSynchronized(false),
	pointNames(),
	pointP(),
	resExt(),
	tangExt(),
	resInt(),
	tangInt(),
	numDofs(0),
	numDofsP(0),
	numDofsI(0),
	numDofsIP(0),
	numDofsB(0),
	localGlob(),
	localGlobI(),
	vOtherInter(),
	kAux(),
	kAuxp(),
	KIB(),
	KBI(),	// UNCOMMENT FOR UNSYMMETRIC MATRICES
	RB(),
	Y0(),
	YB(),
	B0(),
	KIB2(),
	KIB2copia(),
	KBI2(),	// UNCOMMENT FOR UNSYMMETRIC MATRICES
	KBB2(),
	globalLocal(),
	globalLocalI(),
	globalLocalII(),
	localGlobIII(),
	otherLocalGlob(),
	symmFlag(true),
	isKResAlreadyComputed(false)
{}

// ***

Interface::Interface(int a, int b, const int& integratorType, Process* local )
:	lProc(a),
	rProc(b),	
	localProc(local),
	relatedInterfaceNum(0),
	isSynchronized(false),
	pointNames(),
	pointP(),
	resExt(),
	tangExt(),
	resInt(),
	tangInt(),
	numDofs(0),
	numDofsP(0),
	numDofsI(0),
	numDofsIP(0),
	numDofsB(0),
	localGlob(),
	localGlobI(),
	vOtherInter(),
	kAux(),
	kAuxp(),
	KIB(),
	KBI(),		// UNCOMMENT FOR UNSYMMETRIC MATRICES
	RB(),
	Y0(),
	YB(),
	B0(),
	KIB2(),
	KIB2copia(),
	KBI2(),		// UNCOMMENT FOR UNSYMMETRIC MATRICES
	KBB2(),
	globalLocal(),
	globalLocalI(),
	globalLocalII(),
	localGlobIII(),
	otherLocalGlob(),
	symmFlag(true),
	isKResAlreadyComputed(false)
{
	fillPointNames();
	if (integratorType==2)	symmFlag=false;	// Non-symmetric matrices for conserving integration
}

// ***

Interface::~Interface()
{
	/*
	SUPERLU_FREE (ferr);
	SUPERLU_FREE (berr);

	Destroy_Dense_Matrix(&Y0slu);
	*/
}

// ***

void Interface::fillPointNames()
{
	map<string, int, less<string> >::iterator i;
	map<string, Point, less<string> >::iterator j;
	vector<string>::iterator k;

	// Local points
	for(i=localProc->localSharedPoints.begin(); i!= localProc->localSharedPoints.end() ;i++)
	{
		if( (i->second)==rProc)	pointNames.push_back(i->first);
	}
	// Non-local points. They are on the root system of the local process
	for(j=(localProc->rootSys->mPoints).begin(); j!= (localProc->rootSys->mPoints).end() ;j++)
	{
		// Check if the point belongs to a related process
		if( localProc->mSysProc[getStringHead(j->first)]==rProc)	pointNames.push_back(j->first);
	}
	// Remove duplicate entries at vector pointNames
	// First, sort
	sort(pointNames.begin(),pointNames.end());
	// Second, remove consecutive duplicates
	k = unique(pointNames.begin(),pointNames.end());
	pointNames.erase(k,pointNames.end());
	vector<string>(pointNames).swap(pointNames);
//	cout << "at interface " << lProc << "," << rProc << " pointNames es " << endl;
//	printSTLVector(pointNames);
}

// ***

bool Interface::checkRelation(int local, int remote)
{
	if(local==rProc && remote==lProc) return(true); else return(false);
}

// ***

void Interface::setRemoteInterface()
{
	int i;
	int numInter = localProc->getPTProc(rProc)->getNumberInterfaces();
	bool rel;
	for(i = 0; i < numInter ; i++)
	{
		rel = localProc->getPTProc(rProc)->vInterfaces[i].checkRelation(lProc,rProc);
		if(rel==true) break;
	}
	relatedInterfaceNum = i;
}

// ***

void Interface::synchronize()
{
	/** First, check that the related interface is not already synchronized */
	if( ( localProc->getPTProc(rProc)->vInterfaces[relatedInterfaceNum].isSynchronized )==false)
	{
		localProc->getPTProc(rProc)->vInterfaces[relatedInterfaceNum].subsPointNames(pointNames);
		localProc->getPTProc(rProc)->vInterfaces[relatedInterfaceNum].isSynchronized=true;
		isSynchronized=true;
	}
}
// ***

void Interface::subsPointNames(const vector<string>& v)
{
	pointNames = v;
}

// ***

void Interface::getPointers()
{
	pointP = localProc->getPTPoints(pointNames);
}

// ***

void Interface::calcTangResExt(sparseMatrix& Kint, bool b)
{
	int i;
	int j;
	int k;
	double val;
	double val2;

	// Computes tangExt and resExt from KII, KIB, KBI, KBB, RI, RB

	// KII*Q=KIB
	// KII*x=RI
	// Calling Y=( Q | x ) and B=( KIB | RI )
	// Y and B have size numDofsI x (numDofs+1) 
	// the system to be solved is KII*Y=B
	
	if(KIB.size()!=0)	solveKII(Kint, b);	// Provides solution matrices Y0, YB.  Kint is a factorized matrix
	/*
	if (lProc==1)
	{
		cout << "Matriz Y0 en interfaz " << lProc << "," << rProc << " es " << endl;
		for (i=0; i<Y0.size(); ++i)
		{
			for (j=0; j<Y0[i].size(); ++j)	cout << Y0[i][j] << " ";
			cout << endl;
		}
		cout << "Matriz YB en interfaz " << lProc << "," << rProc << " es " << endl;
		for (i=0; i<YB.size(); ++i)
		{
			for (j=0; j<YB[i].size(); ++j)	cout << YB[i][j] << " ";
			cout << endl;
		}		
	}
	*/
	// tangExt = KBI*Q + KBB
	for (i=0; i<numDofs; ++i)
	{
		for (j=0; j<numDofs; ++j)
		{
			val = 0.0;
			for (k=0; k<numDofsIP; ++k)
			{
				if (symmFlag==false)	val2 = KBI[i][ globalLocalII[k] ];	// UNSYMMETRICAL matrix
				else				val2 = -KIB[ globalLocalII[k] ][i];	// SYMMETRICAL matrix
				//** WARNING: Remember that Y0 is used transposed for efficiency
				//** the original line was				if (val2!=0.0) val+= val2*Y0[k][j];
				if (val2!=0.0) val+= val2*Y0[j][k];
			}
			for (k=0; k<numDofsB; ++k)
			{
				if (symmFlag==false)	val2 = KBI[i][ localGlobIII[k] ];	// UNSYMMETRICAL matrix
				else				val2 = -KIB[ localGlobIII[k] ][i];	// SYMMETRICAL matrix
				//** WARNING: Remember that YB is used transposed for efficiency
				//** original line was	if (val2!=0.0) val+= val2*YB[k][j];
				if (val2!=0.0) val+= val2*YB[j][k];
			}
			tangExt[i][j] = val + KBB[i][j];
		}
	}
	/*
	if (lProc==1)
	{
		cout << "Matriz Y0 en interfaz " << lProc << "," << rProc << " es " << endl;
		for (i=0; i<Y0.size(); ++i)
		{
			for (j=0; j<Y0[i].size(); ++j)	cout << Y0[i][j] << " ";
			cout << endl;
		}
		cout << "Matriz YB en interfaz " << lProc << "," << rProc << " es " << endl;
		for (i=0; i<YB.size(); ++i)
		{
			for (j=0; j<YB[i].size(); ++j)	cout << YB[i][j] << " ";
			cout << endl;
		}		
	}
	*/
	
	// resExt = RB-KBI*x
	for (i=0; i<numDofs; ++i)
	{
		val = 0.0;
		for (k=0; k<numDofsIP; ++k)
		{
			if (symmFlag==false)	val2 = KBI[i][ globalLocalII[k] ];	// UNSYMMETRICAL matrix
			else				val2 = -KIB[ globalLocalII[k] ][i];	// SYMMETRIC matrix
			//** WARNING: Remember that Y0 is used transposed for efficiency
			//** the original line was			if (val2!=0.0) val+= val2*Y0[k][numDofs];
			if (val2!=0.0) val+= val2*Y0[numDofs][k];
		}
		for (k=0; k<numDofsB; ++k)
		{
			if (symmFlag==false)	val2 = KBI[i][ localGlobIII[k] ];	// UNSYMMETRIC matrix
			else				val2 = -KIB[ localGlobIII[k] ][i];	// SYMMETRIC matrix			
			//** WARNING: Remember that YB is used transposed for efficiency
			//** original line was	if (val2!=0.0) val+= val2*YB[k][numDofs];
			if (val2!=0.0) val+= val2*YB[numDofs][k];
		}
		resExt[i] = RB[i] - val;
	}
	/*
	if (lProc==1)
	{
		cout << "tangExt is " << tangExt.size() << endl;
		for (i=0; i<tangExt.size(); ++i)
		{
			for (j=0; j<tangExt[i].size(); ++j)	cout << tangExt[i][j] << " ";
			cout << endl;
		}
		cout << "resExt is " << endl;
		printSTLVector(resExt);
	}
	*/
}

// ***

void Interface::initializeKR(vector<vector<double> >& Kbb, vector<vector<double> >& Kib, vector<vector<double> >& Kbi,
								vector<int>& globalLocalInt, vector<int>& globalLocalBou,
								vector<double>& isoRes, sparseMatrix& Kint, int& matType)
{	
	int i;
	int j;
	int k;
	int a;
	int b;
	int c;
	int numDofOther;
	double auxd;

	
	// Initializing KIB, KBI and RI (which is the last column of KIB)
	// just for sparse matrices
	for (i=0; i<numDofsI; ++i)
	{
		a = globalLocalInt[ localGlobI[i] ];
		for (j=0; j<numDofs; ++j)
		{
			b = globalLocalInt[ localGlob[j] ];
			//if (lProc==0) cout << "a = " << a << " , b = " << b << endl;
			if (a==-1 && b==-1)	// Belongs to Kbb
			{
				KIB[i][j] = -Kbb[ globalLocalBou[ localGlobI[i] ] ][ globalLocalBou[ localGlob[j] ] ];
				if (symmFlag==false)	KBI[j][i] = Kbb[ globalLocalBou[ localGlob[j] ] ][ globalLocalBou[ localGlobI[i] ] ];
			}
			else if (a!=-1 && b==-1)	// Belongs to Kib
			if (a!=-1 && b==-1)	// Belongs to Kib
			{
				KIB[i][j] = -Kib[ globalLocalInt[ localGlobI[i] ] ][ globalLocalBou[ localGlob[j] ] ];
				if (symmFlag==false)	KBI[j][i] = Kbi[ globalLocalBou[ localGlob[j] ] ][ globalLocalInt[ localGlobI[i] ] ];
			}
			else if (a==-1 && b!=-1)	// Belongs to Kbi
			{
				if (symmFlag==false)	KIB[i][j] = -Kbi[ globalLocalBou[ localGlobI[i] ] ][ globalLocalInt[ localGlob[j] ] ];
				else					KIB[i][j] = -Kib[ globalLocalInt[ localGlob[j] ] ][  globalLocalBou[ localGlobI[i] ] ];
				
				if (symmFlag==false)	KBI[j][i] = Kib[ globalLocalInt[ localGlob[j] ] ][ globalLocalBou[ localGlobI[i] ] ];
			}
		}
		KIB[i][numDofs] = isoRes[ localGlobI[i] ];
	}
	
	
	
	// Initializing KBB and RB
	// just sparse matrices
	for (i=0; i<numDofs; ++i)
	{
		a = globalLocalInt[ localGlob[i] ];
		RB[i] = isoRes [ localGlob[i] ];
		for (j=0; j<numDofs; ++j)	
		{
			b = globalLocalInt[ localGlob[j] ];
			if (a==-1 && b==-1)	// Belongs to Kbb
			{
				KBB[i][j] = Kbb[ globalLocalBou[ localGlob[i] ] ][ globalLocalBou[ localGlob[j] ] ];
			}
			else if (a!=-1 && b==-1)	// Belongs to Kib
			if (a!=-1 && b==-1)	// Belongs to Kib
			{
				KBB[i][j] = Kib[ globalLocalInt[ localGlob[i] ] ][ globalLocalBou[ localGlob[j] ] ];
			}
			else if (a==-1 && b!=-1)	// Belongs to Kbi
			{
				if (symmFlag==false)	KBB[i][j] = Kbi[ globalLocalBou[ localGlob[i] ] ][ globalLocalInt[ localGlob[j] ] ];
				else					KBB[i][j] = Kib[ globalLocalInt[ localGlob[j] ] ][ globalLocalBou[ localGlob[i] ] ];
			}
		}
	}
	//		cout << "KBB of interface " << lProc << "," << rProc << " already computed " << endl;

	
	// Only if there are other boundaries (numDofsB!=0)
	// Just sparse matrices
	if(numDofsB!=0)
	{	
		// Initializing KIB2
		for (i=0; i<numDofsIP; ++i)
		{
			a = globalLocalInt[ localGlobI[ globalLocalII[i] ] ];
			for (j=0; j<numDofsB; ++j)	
			{
				b = globalLocalInt[ localGlobI[ localGlobIII[j] ] ];
				if (a==-1 && b==-1)	// Belongs to Kbb
				{
					KIB2[j][i] = -Kbb[ globalLocalBou[ localGlobI[ globalLocalII[i] ] ] ][ globalLocalBou[ localGlobI[ localGlobIII[j] ] ] ];
					if (symmFlag==false) KBI2[j][i] = Kbb[ globalLocalBou[ localGlobI[ localGlobIII[j] ] ] ][ globalLocalBou[ localGlobI[ globalLocalII[i] ] ] ];
				}
				else if (a!=-1 && b==-1)
				{
					KIB2[j][i] = -Kib[ globalLocalInt[ localGlobI[ globalLocalII[i] ] ] ][ globalLocalBou[ localGlobI[ localGlobIII[j] ] ] ];
					if (symmFlag==false) KBI2[j][i]  = Kbi[ globalLocalBou[ localGlobI[ localGlobIII[j] ] ] ][ globalLocalInt[ localGlobI[ globalLocalII[i] ] ] ];
				}
				else if (a==-1 && b!=-1)
				{
					if (symmFlag==false)	KIB2[j][i] = -Kbi[ globalLocalBou[ localGlobI[ globalLocalII[i] ] ] ][ globalLocalInt[ localGlobI[ localGlobIII[j] ] ] ];
					else					KIB2[j][i] = -Kib[ globalLocalInt[ localGlobI[ localGlobIII[j] ] ] ][  globalLocalBou[ localGlobI[ globalLocalII[i] ] ] ];
					if (symmFlag==false)	KBI2[j][i] = Kbi[ globalLocalBou[ localGlobI[ localGlobIII[j] ] ] ][ globalLocalInt[ localGlobI[ globalLocalII[i] ] ] ];
				}
			}
		}
		
		//			cout << "Initial KIB2 and KBI2 of interface " << lProc << "," << rProc << " already computed " << endl;
	}
}

// ***

void Interface::setInterExt(vector<vector<double> >& Kbb, vector<vector<double> >& Kib, vector<vector<double> >& Kbi,
							   vector<int>& globalLocalInt, vector<int>& globalLocalBou,
							   vector<double>& isoRes, sparseMatrix& Kint, int& matType, bool bFlag)
{
	if (bFlag==true)	isKResAlreadyComputed=false;	// bFlag=true means that everything needs to be computed, because it is the first Newton iteration 
	
	if (isKResAlreadyComputed==false || vOtherInter.size()>0 )  setInterExtNew(Kbb, Kib, Kbi, globalLocalInt, globalLocalBou, isoRes, Kint, matType, bFlag);
}	
// ***

void Interface::setInterExtNew(vector<vector<double> >& Kbb, vector<vector<double> >& Kib, vector<vector<double> >& Kbi,
								vector<int>& globalLocalInt, vector<int>& globalLocalBou,
								vector<double>& isoRes, sparseMatrix& Kint, int& matType, bool bFlag)
{
	int i;
	int j;
	int k;
	int a;
	int b;
	int c;
	int numDofOther;
	double auxd;
	
	// Initializing KIB, KBI and RI (which is the last column of KIB)
	// just for sparse matrices
	for (i=0; i<numDofsI; ++i)
	{
		KIB[i][numDofs] = isoRes[ localGlobI[i] ];
	}
	
	// KBB and RB already initialized
	// just sparse matrices
	
	// Add to KIB, KBI, RI, KBB and RB internal contributions from other interfaces
	
	//	if (lProc==1)	cout << "A punto de añadir a KIB, KBI, KBB, RI y RB ... " << endl;
	for(k=0; k<vOtherInter.size(); k++)	// Loop on the rest of interfaces
	{
		numDofOther = otherLocalGlob[k].size();	
		for(i=0; i<numDofOther; i++)
		{
			a = globalLocalI [ otherLocalGlob[k][i] ];	// a=-1 means that the global dof is not at the interior
			c = globalLocal [ otherLocalGlob[k][i] ];	// c=-1 means that the global dof is not at the boundary
			if (a!=-1) 
			{
				KIB[a][numDofs] += vOtherInter[k]->resInt[i];
				//				if (lProc==1)	cout << "Added value to the residual RI at interface " << lProc << "," << rProc << endl;
			}
		}
	}
	
	//		cout << "KIB and KBI of interface " << lProc << "," << rProc << " already computed " << endl;
	
	// Only if there are other boundaries (numDofsB!=0)
	// Just sparse matrices
	if(numDofsB!=0)
	{
		//			cout << "Initial KIB2 and KBI2 of interface " << lProc << "," << rProc << " already computed " << endl;
		// Initializing KBB2
		for (i=0; i<numDofsB; ++i)
		{
			a = globalLocalInt[ localGlobI[ localGlobIII[i] ] ];
			for (j=0; j<numDofsB; ++j)
			{
				b = globalLocalInt[ localGlobI[ localGlobIII[j] ] ];
				if (a==-1 && b==-1)	// Belongs to Kbb
				{
					KBB2[i][j] = Kbb[ globalLocalBou[ localGlobI[ localGlobIII[i] ] ] ][ globalLocalBou[ localGlobI[ localGlobIII[j] ] ] ];
				}
			}
		}
		//			cout << "KBB2 of interface " << lProc << "," << rProc << " already computed " << endl;
		// Adding in KIB2 and KBI2 internal contributions from other interfaces	
		
		//		if (lProc==1)	cout << "A punto de añadir a KIB2, KBI2 y KBB2 ... " << endl;
		for(k=0; k<vOtherInter.size(); k++)	// Loop on the rest of interfaces
		{
			numDofOther = otherLocalGlob[k].size();	
			for(i=0; i<numDofOther; i++)
			{
				//	a = localProc->globalLocalI [otherLocalGlob[k][i] ];	// a=-1 means that the global dof is not at the interior of process
				c = globalLocalIII [ globalLocalI [ otherLocalGlob[k][i] ] ];	// c=-1 means that the global dof is not at other boundaries
				for(j=0; j<numDofOther; j++)
				{
					b = globalLocalIII [ globalLocalI [ otherLocalGlob[k][j] ] ];	// b=-1 means that the global dof is not at other boundaries
					if (c!=-1 && b!=-1)
					{
						auxd = vOtherInter[k]->tangInt[i][j];
						if (auxd!=0.0)
						{
							KBB2[c][b]+= auxd;
							//							if (lProc==1)	cout << "Added value to KBB2 at interface " << lProc << "," << rProc << endl;
						}
					}
				}
			}
		}
	}
	
	calcTangResExt(Kint, bFlag);
	
	isKResAlreadyComputed = true;
}

// ***

void Interface::setInterInt()
{
	extrude();
	
	//if (lProc==0)   cout << "tangSize+residual = " << kAux.size() << endl;
	request[0] = MPI::COMM_WORLD.Isend(&kAux[0],tam,MPI::DOUBLE,rProc,0);
	request[1] = MPI::COMM_WORLD.Irecv(&kAuxp[0],tam,MPI::DOUBLE,rProc,0);
	
	request[0].Waitall(2,request,status);
	
	unextrude();
	
	/*
// Non-blocking implementation
	
	request[0] = MPI::COMM_WORLD.Isend(&resExt[0],resExt.size(),MPI::DOUBLE,rProc,0);
	request[1] = MPI::COMM_WORLD.Irecv(&resInt[0],resInt.size(),MPI::DOUBLE,rProc,0);
	
	extrude();

	request[2] = MPI::COMM_WORLD.Isend(&kAux[0],kAux.size(),MPI::DOUBLE,rProc,0);
	request[3] = MPI::COMM_WORLD.Irecv(&kAuxp[0],kAuxp.size(),MPI::DOUBLE,rProc,0);
	
	request[0].Waitall(4,request,status);

	unextrude();
*/
}

// ***

void Interface::resetInterInt()
{
	int i;
	int j;
	for(i=0; i<resInt.size(); i++)
	{
		resInt[i] = 0.0;
		for(j=0; j<resInt.size(); j++)  tangInt[i][j] = 0.0;
	}
}

// ***

void Interface::calcSize()
{
	int i;
	int j;
	
	vector<int> ua;
	
	// Set localGlob
	for(i=0; i<pointP.size(); i++)		// Loop on points
	{
		ua = pointP[i]->globalTNum;		// Contains the global Tdofs numbers of the point 
		for(j=0; j<ua.size(); j++)
		{
			localGlob.push_back(ua[j]) ;
		}
		ua = pointP[i]->globalRNum;		// Contains the global Rdofs numbers of the point 
		for(j=0; j<ua.size(); j++)
		{
			localGlob.push_back(ua[j]) ;
		}
	}
	numDofs = localGlob.size();
	numDofsP = localProc->numDofs;	// Total dofs of the parent process
	numDofsI = numDofsP - numDofs;
}

// ***

void Interface::setSize()
{
	int i;
	int j;
	int k;
	vector<int> auxGlobalNum;	// Auxiliary vector containing the global dof numbers of the interface points
//	vector<int>::iterator it;
	int cont;

    numDofsIP = localProc->numDofsI;	// Total interior dofs of the parent process
	numDofsB = numDofsI - numDofsIP;	// Difference between the interior dofs at interface and process

	// localGlobI. Input: local interior dof number / Output: global (process) dof number
	// First, get the global numbers of the interface dofs, store them at auxGlobalNum
	// Remember that global dof numbers are correlative from 0, stablished for each process
	for( i=0;i<numDofs;i++)
	{
		auxGlobalNum.push_back(localGlob[i]);
	}	
	for( i=0; i<numDofsP; ++i)
	{
		// If the global dof number is not at interface, it is at the interior
		if ( count (auxGlobalNum.begin(), auxGlobalNum.end(), i)==0 ) localGlobI.push_back(i);
	}

	// globalLocalI.  Input: global (process) dof number / Output: local interior dof number
	// If the global number is not in the interior, put -1
	// Perform this search with stl. TO BE DONE.
	globalLocalI.resize(numDofsP);
	for (i=0; i<numDofsP; ++i)
	{
		cont = 0;
		for (j=0; j<numDofsI; ++j)
		{
			if ( i== localGlobI[j] )
			{
				cont++;
				break;
			}
		}
		if (cont!=0) globalLocalI[i]=j; else globalLocalI[i]=-1;
	}
	
	// globalLocalII. Input: global (process) dof number / Output: local interior_interior dof number
	globalLocalII.resize(numDofsIP);
	for (i=0; i<numDofsIP; ++i)
	{
		globalLocalII[i] = globalLocalI[ localProc->localGlobI[i] ];
	}
	
	// localGlobIII. Input: local other-boundaries dof / Output: interior dof in interface
	for (i=0; i<numDofsI; ++i)
	{
		cont=0;
		for (j=0; j<numDofsIP; ++j)
		{
			if ( globalLocalII[j] == i )
			{
				cont++;
				break;
			}
		}
		if (cont==0) localGlobIII.push_back(i);
	}

	// globalLocalIII. Input: interior dof in interface / Output: local other-boundaries dof
	globalLocalIII.resize(numDofsI);
	for (i=0; i<numDofsI; ++i)
	{
		cont = 0;
		for (j=0; j<numDofsB; ++j)
		{
			if ( i== localGlobIII[j] )
			{
				cont++;
				break;
			}
		}
		if (cont!=0) globalLocalIII[i]=j; else globalLocalIII[i]=-1;
	}
	
	// globalLocal. Input: Global dof at process / Output: local boundaries dof
	globalLocal.resize(numDofsP);
	for (i=0; i<numDofsP; ++i)
	{
		cont = 0;
		for (j=0; j<numDofs; ++j)
		{
			if ( i== localGlob[j] )
			{
				cont++;
				break;
			}
		}
		if (cont!=0) globalLocal[i]=j; else globalLocal[i]=-1;
	}

	// otherLocalGlob.  Each row is the localGlob vector of the other interfaces of the process
	for(k=0; k<vOtherInter.size(); k++)	// Loop on the rest of interfaces
	{
		otherLocalGlob.push_back (vOtherInter[k]->localGlob);
	}
		
	resExt.resize(numDofs);		// Resize external res
	tangExt.resize(numDofs);	// Resize external tang
	for( i=0;i<numDofs;i++) tangExt[i].resize(numDofs);

	resInt.resize(numDofs);		// Resize internal res
	tangInt.resize(numDofs);	// Resize internal tang
	for( i=0;i<numDofs;i++) tangInt[i].resize(numDofs);
	
//	tam = (int)	numDofs*(1+numDofs)/2 ;
	tam = numDofs*numDofs + numDofs;
	kAux.resize( tam); 	// Resize auxiliary vector used in extrude, for MPI communications
	kAuxp.resize( tam); 	// Resize auxiliary vector used in extrude, for MPI communications

	KIB.resize(numDofsI);	// Tangent of interior - boundary dofs. Last column is for RI
	for( i=0;i<numDofsI;i++) KIB[i].resize(numDofs+1);

	if (symmFlag==false)
	{
		KBI.resize(numDofs);	// Tangent of boundary - interior dofs UNCOMMENT FOR UNSYMMETRICAL MATRICES
		for( i=0;i<numDofs;i++) KBI[i].resize(numDofsI);
	}

	KBB.resize(numDofs);	// Tangent of boundary - boundary dofs
	for( i=0;i<numDofs;i++) KBB[i].resize(numDofs);
	
	RB.resize(numDofs);
	
//	cout << "En interface " << lProc << " , " << rProc << " numDofsI es " << numDofsI << " y numDofs es " << numDofs << endl;
	
	// WARNING: Y0 is the transposed of the actual matrix, for efficiency when solving with Kint
	// Original lines were:
	//*	Y0.resize(numDofsIP);
	//*	for( i=0;i<numDofsIP;i++) Y0[i].resize(numDofs+1);

	Y0.resize(numDofs+1);
	for( i=0;i<numDofs+1;i++) Y0[i].resize(numDofsIP);
	Y0ori.resize(numDofs+1);
	for( i=0;i<numDofs+1;i++) Y0ori[i].resize(numDofsIP);


	
	// WARNING: B0 is the transposed of the actual matrix, for efficiency when solving with Kint
	// Original lines were:
	//*	B0.resize(numDofsIP);
	//*	for(i=0; i<numDofsIP;i++)	B0[i].resize(numDofs+1);

	B0.resize(numDofs+1);
	for(i=0; i<numDofs+1;i++)	B0[i].resize(numDofsIP);

	if (numDofsB!=0)
	{
		
		//** WARNING: YB is the transposed of the actual matrix, for efficiency 
		//** Original lines were:
		/*
		YB.resize(numDofsB);
		for(i=0; i<numDofsB;i++)	YB[i].resize(numDofs+1);
		*/
		YB.resize(numDofs+1);
		for(i=0; i<numDofs+1;i++)	YB[i].resize(numDofsB);
//		YBori.resize(numDofs+1);
//		for(i=0; i<numDofs+1;i++)	YBori[i].resize(numDofsB);

		// WARNING: KIB2, KIB2copia are the transposed of the actual matrices, for efficiency when solving with Kint
		// Original lines were:
		/*			//*
		 KIB2.resize(numDofsIP);
		 for(i=0; i<numDofsIP;i++)	KIB2[i].resize(numDofsB);
		 
		 KIB2copia.resize(numDofsIP);
		 for(i=0; i<numDofsIP;i++)	KIB2copia[i].resize(numDofsB);
		 */		
		KIB2.resize(numDofsB);
		for(i=0; i<numDofsB;i++)	KIB2[i].resize(numDofsIP);
		
		KIB2copia.resize(numDofsB);
		for(i=0; i<numDofsB;i++)	KIB2copia[i].resize(numDofsIP);		
		
		if (symmFlag==false)
		{
			KBI2.resize(numDofsB);	// UNCOMMENT FOR UNSYMMETRICAL MATRICES
			for(i=0; i<numDofsB;i++)	KBI2[i].resize(numDofsIP);	// WARNING: It is not transposed
		}

		KBB2.resize(numDofsB);
		for(i=0; i<numDofsB;i++)	KBB2[i].resize(numDofsB);
//		KBB2copia.resize(numDofsB);
//		for(i=0; i<numDofsB;i++)	KBB2copia[i].resize(numDofsB);
	}
/*
	cout << "globalLocal en interface " << lProc << "," << rProc << endl;
	printSTLVector(globalLocal);
	cout << "globalLocalI en interface " << lProc << "," << rProc << endl;
	printSTLVector(globalLocalI);
	cout << "globalLocalII en interface " << lProc << "," << rProc << endl;
	printSTLVector(globalLocalII);
	cout << "globalLocalIII en interface " << lProc << "," << rProc << endl;
	printSTLVector(globalLocalIII);
	cout << "localGlobIII en interface " << lProc << "," << rProc << endl;
	printSTLVector(localGlobIII);
	
	cout << "otherLocalGlob en interface " << lProc << "," << rProc << endl;
	printSTLMatrix(otherLocalGlob);
	cout << "KBBlocate en interface " << lProc << "," << rProc << endl;
	printSTLMatrix(KBBlocate);
	cout << "KBB2locate en interface " << lProc << "," << rProc << endl;
	printSTLMatrix(KBB2locate);
	cout << "KIBlocate en interface " << lProc << "," << rProc << endl;
	printSTLMatrix(KIBlocate);
	cout << "KIB2locate en interface " << lProc << "," << rProc << endl;
	printSTLMatrix(KIB2locate);
	cout << "KBIlocate en interface " << lProc << "," << rProc << endl;
	printSTLMatrix(KBIlocate);
	cout << "KBI2locate en interface " << lProc << "," << rProc << endl;
	printSTLMatrix(KBI2locate);
*/
/*
	int n = numDofsIP;
	int nrhs = numDofs+1;
	
	if ( !(rhsx = doubleMalloc(n * nrhs)) ) ABORT("Malloc fails for rhsx[].");
	if ( !(rhsb = doubleMalloc(n * nrhs )) ) ABORT("Malloc fails for rhsb[].");
*/
	/*
	if ( !(ferr = (double *) SUPERLU_MALLOC(nrhs * sizeof(double))) )	
		ABORT("SUPERLU_MALLOC fails for ferr[].");
	
	if ( !(berr = (double *) SUPERLU_MALLOC(nrhs * sizeof(double))) )
		ABORT("SUPERLU_MALLOC fails for berr[].");
	dCreate_Dense_Matrix(&Y0slu, n, nrhs, rhsb, n, SLU_DN, SLU_D, SLU_GE);
//	dCreate_Dense_Matrix(&X, n, nrhs, rhsx, n, SLU_DN, SLU_D, SLU_GE);
	*/
}

//

void Interface::addRes(vector<double>& res_in)
{
	int i;
	double auxd;
	for(i=0; i<numDofs; i++)
	{
		auxd = resInt[i];
		if(auxd!=0.0) res_in[ localGlob[i] ] = res_in[ localGlob[i] ] + auxd;
	}
}

// ***

void Interface::addTang(vector<vector<double> >& Kbb, vector<vector<double> >& Kib, vector<vector<double> >& Kbi,
						vector<int>& globalLocalInt, vector<int>& globalLocalBou,
						int& matType)
{
	// Just for sparse matrices
	
	int i;
	int j;
	double auxd;
//	int a, b;
		
	for(i=0; i<numDofs; i++)
	{
//		a = globalLocalInt[ localGlob[i] ];
		for(j=0; j<numDofs; j++)
		{
//			b = globalLocalInt[ localGlob[j] ];
			auxd =  tangInt[i][j];
			if (auxd!=0.0)
			{
//				if (a==-1 && b==-1)	// Belongs to Kbb
				{
					Kbb[ globalLocalBou[ localGlob[i] ] ][ globalLocalBou[ localGlob[j] ] ]+= auxd;
//					if (lProc==1)	cout << "Added value to Kbb at interface " << lProc << "," << rProc << endl;
				}
				/*
				else if (a!=-1 && b==-1)	// Belongs to Kib
				{
					Kib[ globalLocalInt[ localGlob[i] ] ][ globalLocalBou[ localGlob[j] ] ]+= auxd;
//					if (lProc==1)	cout << "Added value to Kib at interface " << lProc << "," << rProc << endl;
				}
				else if (a==-1 && b!=-1)	// Belongs to Kbi
				{
					if(symmFlag==false)	Kbi[ globalLocalBou[ localGlob[i] ] ][ globalLocalInt[ localGlob[j] ] ]+= auxd;
				}
				*/
			}
		}
	}
}

// ***

void Interface::setOtherInter(vector<Interface*> v)
{
	vOtherInter = v;
}

// ***

void Interface::extrude()
{
	int i;
	int j;

	if (symmFlag==false)
	{
		// Without assuming symmetry
		for(i=0; i<numDofs; i++)
		{
			for(j=0; j<numDofs; j++)
			{
				kAux[numDofs*i+j] = tangExt[i][j];
			}
		}
	}
	else
	{
		// Assuming symmetry
		int cont = 0;
		for (i=0; i<numDofs; ++i)
		{
			for (j=i; j<numDofs; ++j) 
			{
				kAux[cont] = tangExt[i][j];
				cont+=1;
			}
		}
	}
	for (i=0; i<numDofs; ++i)
	{
		kAux[numDofs*numDofs+i] = resExt[i];
	}	
}

// ***

void Interface::unextrude()
{
	int i;
	int j;
	
	if (symmFlag==false)
	{
		// Without assuming symmetry
		for(i=0; i<numDofs; i++)
		{
			for(j=0; j<numDofs; j++)
			{
				tangInt[i][j] = kAuxp[numDofs*i+j];
			}
		}
	}
	else
	{
		// Assuming symmetry
		int cont = 0;
	
		for (i=0; i<numDofs; ++i)
		{
			for (j=i; j<numDofs; ++j) 
			{
				tangInt[i][j] = kAuxp[cont];
				cont+=1;
			}
		}
		for (i=0; i<numDofs; ++i)
		{
			for (j=0; j<i; ++j) 
			{
				tangInt[i][j] = tangInt[j][i];
			}
		}		
	}
	for (i=0; i<numDofs; ++i)
	{
		resInt[i] = kAuxp[i+numDofs*numDofs];
	}	
}

// ***

bool Interface::checkIfDof(const int num)
{
	if ( count(localGlob.begin(), localGlob.end(), num) == 0 ) return false; else return true;
}

// ***

void Interface::solveKII(sparseMatrix& Kint, bool b)
{	
	// The input matrix Kint is already factorized
	
//	cout << "In solveKII of interface " << lProc << "," << rProc << " ...\n";
	
	int i;
	int j;
	int k;
	double val;
	double val2;

//	cout << "b flag is " << b << endl;
	if(b==true)
	{
		// B0
		setZero(B0);

		for (i=0; i< numDofsIP; ++i)
		{
			for (j=0; j<numDofs+1; ++j)
			{
				val = KIB[ globalLocalII[i] ][j];
				//** WARNING: remenber that B0 is used transposed
				//** original line was			if (val!=0.0) B0[i][j] = val;
				if (val!=0.0) B0[j][i] = val;
			}	
		}

	
		Y0ori = B0;
		// For openMP implementation
		/*
		if (Kint.rows()!=0)
		{
			if (!Kint.typeName.compare("slu"))
			{
				B_temp = convert_stl_SLU(B0);
				Kint.subsSolveTrans(B_temp, B_temp);	// Solution of Kint*X=B0 is in B0
				convert_SLU_stl(B_temp,B0);
				Destroy_Dense_Matrix(&B_temp);
			}
			else 
			{
				Kint.subsSolveTrans(B0);	// Solution of Kint*X=B0 is in B0
			}
		}
		*/
		if (Kint.rows()!=0)
		{
			Kint.subsSolveTrans(B0);	// Solution of Kint*X=B0 is in B0
		} 
	}
	/*
	if (lProc==1)
	{
		cout << "X dentro de interface DESPUES" << lProc << "," << rProc << " es" << endl;
		for (i=0; i<B0.size(); ++i)
		{
			for (j=0; j<B0[i].size(); ++j)	cout << B0[i][j] << " ";
			cout << endl;
		}
	}
	*/
	// Only if there are other boundaries (numDofsB!=0)
	if (numDofsB!=0)
	{
		// BB.  Stored at the begining at YB
		setZero(YB);
		for (i=0; i<numDofsB; ++i)
		{
			for (j=0; j<numDofs+1; ++j)
//			if (b==true)
//			{
//				for (j=0; j<numDofs; ++j)
				{
					val = KIB[ localGlobIII[i] ][j];
					//** WARNING: remenber that YB is used transposed
					//** original line was	if (val!=0.0) YB[i][j] = val;
					if (val!=0.0) YB[j][i] = val;
				}
//			}
//			val = KIB[ localGlobIII[i] ][numDofs];
//			if (val!=0.0) YB[numDofs][i] = val; else YB[numDofs][i]=0.0;
		}

		if (b==true)
		{
			KIB2copia = KIB2;
			// For openMP implementation
			/*
			if (Kint.rows()!=0)
			{
				if (!Kint.typeName.compare("slu"))
				{
					B_temp = convert_stl_SLU(KIB2copia);
					Kint.subsSolveTrans(B_temp, B_temp);	// Solution of Kint*Q=-KIB2 is in KIB2copia=Q
					convert_SLU_stl(B_temp,KIB2copia);
					Destroy_Dense_Matrix(&B_temp);
				}
				else
				{
					Kint.subsSolveTrans(KIB2copia);	// Solution of Kint*Q=-KIB2 is in KIB2copia=Q
				}				
			}
			*/
			if (Kint.rows()!=0)
			{
				Kint.subsSolveTrans(KIB2copia);	// Solution of Kint*Q=-KIB2 is in KIB2copia=Q
			}
		}
		//JC KIB2copia = KIB2;
			
		//JC if (Kint.rows()!=0)
		/*
		if (Kint.rows()!=0 && b==true)
		{
			Kint.subsSolveTrans(KIB2copia);	// Solution of Kint*Q=-KIB2 is in KIB2copia=Q
		}
		*/
/*		
		cout << "Q' dentro de interface " << lProc << "," << rProc << " es" << endl;
		printSTLMatrix(KIB2);

		cout << "KBB2 ANTES dentro de interface " << lProc << "," << rProc << " es" << endl;
		printSTLMatrix(KBB2);
*/
		//JC
		//KBB2copia = KBB2;
		
		// k = KBB2 + KBI2*Q
		for (i=0; i<numDofsB; ++i)
		{
			for (j=0; j<numDofsB; ++j)
			{
				val = 0.0;
				for (k=0; k<numDofsIP; ++k)
				{
					//JC if (symmFlag==false)	val+= KBI2[i][k]*KIB2[j][k];		// UNSYMMETRIC matrix
					//JC else					val+= -KIB2copia[i][k]*KIB2[j][k];	// SYMMETRIC matrix
					if (symmFlag==false)	val+= KBI2[i][k]*KIB2copia[j][k];		// UNSYMMETRIC matrix
					else					val+= -KIB2[i][k]*KIB2copia[j][k];	// SYMMETRIC matrix			
					//**WARNING: remember that KIB2 is used transposed for efficiency
				}
				if(val!=0.0) KBB2[i][j]+= val;
				//JCif(val!=0.0) KBB2copia[i][j]+= val;
			}
		}
//		cout << "KBB2 dentro de interface " << lProc << "," << rProc << " es" << endl;
//		printSTLMatrix(KBB2);

		
		// r = BB - KBI2*X      r=YB at this moment
		// YB = BB;
		for (i=0; i<numDofsB; ++i)
		{
			for (j=0; j<numDofs+1; ++j)
			{
				val = 0.0;
				for (k=0; k<numDofsIP; ++k)
				{
					//**WARNING: remember that KIB2copia is used transposed for efficiency
					//** original line was				val+= -KIB2copia[k][i]*B0[j][k];
					//JC if (symmFlag==false)	val+= KBI2[i][k]*B0[j][k];			// UNSYMMETRICAL matrix
					//JC else					val+= -KIB2copia[i][k]*B0[j][k];	// SYMMETRICAL matrix
					if (symmFlag==false)	val+= KBI2[i][k]*B0[j][k];			// UNSYMMETRICAL matrix
					else					val+= -KIB2[i][k]*B0[j][k];	// SYMMETRICAL matrix
				}
				//** WARNING: remember that YB is used transposed
				//** original line was	YB[i][j]-= val;
				YB[j][i]-= val;
			}
		}
		
		// kmat*YB=r.
/*
		cout << "k' dentro de interface " << lProc << "," << rProc << " es" << endl;
		printSTLMatrix(KBB2);
		cout << "r' dentro de interface " << lProc << "," << rProc << " es" << endl;
		printSTLMatrix(YB);
*/
		lu(KBB2);	// THIS THE ONLY FACTORIZATION PERFORMED AT EVERY INTERFACE ITERATION
		forSubTrans(KBB2,YB);
		backSubTrans(KBB2,YB);
	
		/*JC
		lu(KBB2copia);	// THIS THE ONLY FACTORIZATION PERFORMED AT EVERY INTERFACE ITERATION
		forSubTrans(KBB2copia,YB);
		backSubTrans(KBB2copia,YB);		
		*/		

	}	// End of the case numDofsB!=0
	
	for (i=0; i<numDofsIP; ++i)
	{
		for (j=0; j<numDofs+1; ++j)
		{
			val = 0.0;
			for (k=0; k<numDofsB; ++k)
			{
				//**WARNING: remember that KIB2 is used transposed for efficiency
				//** original line was				val2 = KIB2copia[i][k];
				//JC val2 = KIB2copia[k][i];
				val2 = KIB2[k][i];
				//** WARNING: remenber that YB is used transposed
				//** original line was	if (val2!=0) val+= val2 * YB[k][j];
				if (val2!=0) val+= val2 * YB[j][k];
			}
			//**WARNING: remember that Y0 is used transposed for efficiency
			//** original line was			Y0[i][j] += val;
//			Y0[j][i] += val;
			Y0[j][i] = Y0ori[j][i] + val;
		}
	}
	/*
	if (Kint.rows()!=0)
	{
		Kint.subsSolveTrans(Y0);	// Solution of Kint*Y0=B0-KIB2*YB
	}
	*/
	
	// Fill right hand side vector (from a matrix)
	
	double t1 = MPI::Wtime();
	// Solution of Kint*Y0=B0-KIB2*YB
	// Solution of Kint*Y01 = B0 is already in B0
	if (KIB2.size()!=0)
	{
		// for openMP implementation
		/*
		if (Kint.rows()!=0)
		{
			if (!Kint.typeName.compare("slu"))
			{
				B_temp = convert_stl_SLU(Y0);
				Kint.subsSolveTrans(B_temp, B_temp);	// Solution of Kint*X=B0 is in B0
				convert_SLU_stl(B_temp,Y0);
				Destroy_Dense_Matrix(&B_temp);
			}
			else
			{
				Kint.subsSolveTrans(Y0);
			}
		}
		*/
		if (Kint.rows()!=0)
		{
			Kint.subsSolveTrans(Y0);
		}
	}
	else Y0 = B0;
	double t2 = MPI::Wtime();	
	
//	cout << "time for computing Y0 at interface " << lProc << "," << rProc << " is " << t2-t1 <<  endl;	
//	cout << "solveKII of interface " << lProc << "," << rProc << " done \n";
}

// ***

void Interface::setSizeExp()
{
	resExt.resize(numDofs);		// Resize external res
	
	resInt.resize(numDofs);		// Resize internal res	
	
//	cout << "Dando tamaño a interface " << lProc << " , " << rProc << endl;
}

// ***

void Interface::initializeRExp(vector<double>& isoRes)
{
	int i;
	for (i=0; i<numDofs; ++i)
	{
		/*
		 if (lProc==1)	
		 {
		 cout << "isoRes.size() = " << isoRes.size() << endl;
		 cout << "localGlobI.size " << localGlob.size() << endl;
		 }
		 */
		resExt[i] = isoRes [ localGlob[i] ];
	}
	
	/*
	 cout << "resExt en proceso " << lProc << endl;
	 for (i=0; i<resExt.size(); ++i)	cout << "  " << resExt[i];
	 cout << endl;
	 */
}

// ***

void Interface::setInterIntExp()
{
	
	// Non-blocking implementation
	
	request[0] = MPI::COMM_WORLD.Isend(&resExt[0],resExt.size(),MPI::DOUBLE,rProc,0);
	request[1] = MPI::COMM_WORLD.Irecv(&resInt[0],resInt.size(),MPI::DOUBLE,rProc,0);
	
	request[0].Waitall(2,request,status);
}

// ***

void Interface::addResExp(vector<double>& res_in)
{
	int i;
	double auxd;
	for(i=0; i<numDofs; i++)
	{
		auxd = resInt[i];
		if(auxd!=0.0) res_in[ localGlob[i] ] = res_in[ localGlob[i] ] + auxd;
	}
}

// ***

void Interface::synchronizeMass(vector<double>& mLumped)
{
	/*
	Process *pRProc = localProc->getPTProc(rProc);
	vector<double> remoteVLumped = pRProc->vInterfaces[relatedInterfaceNum].mLumpedLocal;

	int i;
	
	for (i=0; i<numDofs; ++i)
	{
		mLumped[ localGlob[i] ] += remoteVLumped[i];
	}
	*/

	Process *pRProc;
	//vector<double> remoteVLumped;
	vector<int> remoteLocalGlob;
	
	int i;
	
	// First, check that the related interface is not already synchronized
	if( ( localProc->getPTProc(rProc)->vInterfaces[relatedInterfaceNum].isMassSynchronized )==false)
	{
		//		cout << "Synchronizing interface " << lProc << " , " << rProc << endl;
		pRProc = localProc->getPTProc(rProc);
		//remoteVLumped = pRProc->vLumpedMass;
		remoteLocalGlob = localProc->getPTProc(rProc)->vInterfaces[relatedInterfaceNum].localGlob;
		for (i=0; i<numDofs; ++i)
		{
			//mLumped[ localGlob[i] ] += remoteVLumped[ remoteLocalGlob[i] ];
			mLumped[ localGlob[i] ] += pRProc->vLumpedMass[ remoteLocalGlob[i] ];
			pRProc->vLumpedMass[ remoteLocalGlob[i] ] = mLumped[ localGlob[i] ];
		}
		isMassSynchronized=true;
		localProc->getPTProc(rProc)->vInterfaces[relatedInterfaceNum].isMassSynchronized = true;
	}
}